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ABSTRACT 

The stability of classical semi-implicit scheme, and some more advanced iterative 
schemes recently proposed for NWP purpose is examined. In all these schemes, the 
solution of the centred-implicit non-linear equation is approached by an iterative fixed- 
point algorithm, preconditioned by a simple, constant in time, linear operator. A general 
methodology for assessing analytically the stability of these schemes on canonical prob- 
lems for a vertically unbounded atmosphere is presented. The proposed method is valid 
for all the equation systems usually employed in NWP. However, as in earlier studies, the 
method can be applied only in simplified meteorological contexts, thus overestimating the 
actual stability that would occur in more realistic meteorological contexts. The analy- 
sis is performed in the spatially-continuous framework, hence allowing to eliminate the 
spatial-discretisation or the boundary conditions as possible causes of the fundamental in- 
stabilities linked to the time-scheme itself. The general method is then shown concretely 
to apply to various time-discretisation schemes and equation-systems (namely shallow- 
water, and fully compressible Euler equations). Analytical results found in the literature 
are recovered from the proposed method, and some original results are presented. 
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1 Introduction 

The classical semi-implicit (SI) technique (Robert et ai., 1972) has been widely used 

in NWP since it provides efficient and simple algorithms, at least for spectral models. 

This classical SI method requires the definition of a constant in time linear "reference" 

operator £*, which usually consists in the linearisation of the original system AI, around a 

stationary reference-state, noted X*. For a given state X of the atmosphere, the evolution 

of the system, (dX/dt) = A4.X, is then time-discretised through: 

^ = {M-C*).X + C*:\xf (1) 
ot 

where {6/6t) is the discretised time-derivative operator, and [ ] is the implicit-centred 
temporal average operator. The terms linked to the reference operator C* are thus treated 
in a centred-implicit way, whilst the residual non-linear terms are treated explicitly. For 
this scheme, there is no formal proof of the stability in real atmospheric conditions, due 
to the explicit treatment of non-linear residuals. This prompted the authors of pioneering 
NWP applications of the SI scheme to examine theoretically its stability in idealised 
contexts. 

In a seminal study following this approach, Simmons et al., 1978 (SHB78 hereafter) 
analysed the stability of the SI scheme for the hydrostatic primitive equations (HPE) 
system with a Leap- Frog (3-TL hereafter) time-discretisation by considering the linearised 
equations around a stationary state X (referred to as "atmospheric state" hereafter) when 
the resulting linear "atmospheric" operator £, deviates from the linear "reference" operator 
£* of the SI scheme, thus generating potentially unstable explicitly-treated residuals. 

In the vertically-continuous context they performed a stability analysis valid when the 
eigenfunctions of £ and £* are identical. They showed that when the atmospheric and 
reference temperature profiles (respectively T and T*) are isothermal, the stability of the 
SI scheme requires: 
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< T < 2T*, 



(2) 



hence T* cannot be chosen arbitrarily for applying the 3-TL SI scheme to the HPE system. 

In the finite-difference vertically-discretised context, thew showed that a "vertically- 
discretised analysis" of stability following the same principle simply resulted in the solution 
of a standard eigenvalue problem. They found empirically that a large static-stability for 
the reference-state is necessary to maintain the stability of the scheme for realistic thermal 
atmospheric profiles. As a consequence, they recommended to use a warm isothermal state 
as reference-state, a rule which was then widely adopted for NWP applications using SI 
schemes. 

Finally, they examined the effect of applying a second-order time-filter in the temporal 
average of linear terms, and found that an improved stability is obtained, but at the 
expense of an increased misrepresentation of the wave propagation. 

Cote et al., 1983 (CBS83 hereafter), still in the HPE context, examined the stability of 
the 3-TL SI scheme for a finite-element vertical discretisation using the above vertically- 
discretised analysis method. They established a stability criterion for the 3-TL SI scheme 
in terms of the atmospheric and reference static stabilities (7 and 7* respectively): 



therefore generalizing Q to not necessarily isothermal thermal profiles. 

Still in the HPE context, Simmons and Temperton, 1997 showed with the same method 
that extrapolating two-time level (2-TL) schemes have more stringent stability constraints 
than their 3-TL counterpart. For instance, in the isothermal framework of the SHB78 
analysis, the stability of the 2-TL SI scheme requires: 



< 7 < 27*, 



(3) 



< T < T*. 



(4) 
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As a consequence, they recommended to use a warmer reference temperature than in 
the 3-TL case. ST97 also showed that the 2-TL SI scheme was intrinsically damping 
when stable, a characteristic which was not present in the 3-TL SI scheme. However, 
3-TL schemes require a time-filter to damp the unstable computational modes, and this 
time-filter also damps the transient physical modes. They argue that all these aspects 
being considered, the efi^ective damping of 2-TL and 3-TL schemes is of comparable overall 
intensity. This particular debate will be ignored in this paper, in order to focus on stability 
aspects only. 

The relevance of the SI method for solving numerically the fully compressible Euler 
Equations (EE) was then advocated (Tanguay et al., 1990), and some numerical models in 
EE using this technique were effectively developed: Caya and Laprise (1999), presented a 
model in EE with a 3-TL SI scheme (with a moderate time-decentering first-order accurate 
in time). Semazzi et al. (1995) and Qian et al. (1998) also showed a model in EE but 
with a 2-TL SI scheme (with a strong time-filter however). Besides, the need of more 
robust schemes than the classical SI one for solving the EE system became recognized 
(e.g. Cote et al. 1998, BHBG95), probably motivated by some pathological behaviours 
with the classical SI scheme under some circumstances. 

Schemes with evolution terms treated in a more centred-implicit way are usually be- 
lieved to have an increased robustness, hence fulfilling the latter emerging need, and some 
of them were developed for fine-scale models in the EE system. Bubnova et al., 1995 
used a 3-TL scheme in which the leading non-linear terms of the EE system are treated 
in a centred-implicit way, through a partially iterative method. For a 2-TL scheme. Cote 
et al, 1998, used a fully iterative method, aiming to treat all the evolution terms of the 
EE system in a centred-implicit way. Cullen (2000) examined the benefit of using such 

a fully iterative scheme for the HPE system, arguing that an improved accuracy could 
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be obtained beside the improved stability (this latter being not strictly required however 
for current HPE applications). As a formal justification, he examined the stability of 
this iterative scheme for the 2-TL shallow-water (SW) system. The analysis was limited 
to a scheme called "predictor/corrector", which consists in a single additional iteration 
after the SI scheme. The salient result was that the additional iteration in the "predic- 
tor/corrector" scheme allows to recover a extended range of stability as in Q instead of 
(IH) . In the following, these fully iterative schemes with a more centred-implicit treatment 
of the evolution terms will be referred to as "iterative centred-implicit" (ICI) schemes. 

From the theoretical point of view, the current situation is that no stability analysis 
has been provided for the EE system with SI scheme, and for ICI schemes with more 
iterations, stability analyses are available only for the SW system. 

Here we present a general method to carry out space-continuous stability analyses 
of the various time-discretisation schemes mentioned above, for any usual meteorological 
system of NWP interest (SW, HPE, EE), on canonical problems similar to those examined 
in SHB78. Some original results concerning the EE system and iterative schemes are 
presented. 

This work may also be viewed as a first theoretical investigation about the suitability 
of various time-discretisation schemes for solving numerically the EE system. 

2 General framework for analyses 

The general framework for the stability analyses presented here is basically the same as 
in most earlier studies: The flow is assumed adiabatic inviscid and frictionless in a non- 
rotating dry perfect-gas atmosphere with a Cartesian coordinate system. Moreover, the 
flow is assumed linear around an "atmospheric" basic-state X. The actual evolution of 

the atmospheric flow is thus described by C, the linear-tangent operator to around X. 
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The atmospheric basic-state X is chosen stationary, resting, horizontally homogeneous, 
and hydrostatically balanced. The governing equation for the flow is then: 

where X' = X — X , and the primes are dropped henceforth for clarity. 

Following the usual practice in NWP, the linear operator £* in Jl]) is taken as the 
tangent-linear operator to M. around a reference-state X* which is also chosen stationary, 
resting, horizontally homogeneous, and hydrostatically balanced. Since A" is a resting 
state, the linear Lagrangian time-derivative coincides with the Eulerian time-derivative, 
and the LHS operator of ^ hence holds for Eulerian models as well as for semi-Lagrangian 
models . 



3 The class of ICI schemes in the Hnear framework 

In the restricted resting and linear framework of section [2l classical SI schemes as well 
as the iterative schemes mentioned in the introduction can be gathered in a single class 
of ICI schemes, differing only by their number of iterations. These ICI scheme are first 
presented for a 2-TL discretisation. In this case, the fully implicit-centred (FIC) scheme 
writes: 

X^ - X^ X.X+ +Z.X^ 

-^r- = — 2 — 

where, according to a standard practice for time-discretised equations, the superscripts 
" + " and "0" stand for time levels (t + At) and t respectively. The principle of ICI 
schemes is to approach the FIC solution by starting from an initial "guess" noted X'^^'^\ 
then iterating the following algorithm: 



At _ 2 _ 2 

2 ^ 2 
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(7) 
(8) 



for n = 1,2, Niter- The state, valid a {t + At) is then taken as the last iterated state 
;Y+{Nitcv)_ An examination of this scheme for a model with a single prognostic variable 
without spatial dependency shows that it acts as a fixed-point algorithm for solving the 
implicit non-linear scalar equation f{x) = x, by using an estimate /'* of the derivative /' 
as a preconditioner for convergence. The method converges if |(/' — f'*)/f'*\ < 1. This 
is a weaker condition than the one for the classical (i.e. not preconditioned) fixed-point 
method: |/'| < 1. The initial guess X^^^^ is arbitrary in ICI schemes, but choosing an 
appropriate initial guess may help decreasing the magnitude of the discrepancy between 
FIC and ICI schemes after a fixed number of iterations. 
For 3-TL schemes, the ICI scheme can be defined by: 



2At _ 2 _ 2 



(9) 
(10) 

where the superscript "-" denotes a variable taken at the time-level (t — At). 

Here follows a list of some schemes proposed in the literature, with the corresponding 
characteristics {X'^^'^^ ^Nuer), in the restricted framework of section O 

- Classical 2-TL SI extrapolating scheme: A^^iter = 1 and = {2X° - X~). 

- 2-TL non-extrapolating SI scheme: Nucr = 1 and X'^^^^ = X^. However, this scheme 
is not used in practice since it is only first-order accurate in time, as mentioned in 
Cullen (2000). 

- "Predictor/corrector" scheme of Cullen (2000): iVner = 2 and = X^. 

- Iterative scheme of Cote et al., 1998: general iterative ICI scheme, but used with 

iViter = 2 in practice. The choice of X^^'^^ is not explicitly indicated. 
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- FIC scheme: iViter = oo (does not depend on the choice of X^^^'). This scheme can 
not be achieved in practice for numerical models, but it may be useful for theoretical 
examination of the asymptotic behaviour of the ICI schemes. 

For 3-TL discretisations, the SI scheme, which corresponds to an ICI scheme with Alitor = 1 
and = {2X^ — X~) is the only one to be used in practice. However, iterated 3-TL 

ICI schemes could be used as well, and the 3-TL FIC scheme is equivalent to a 2-TL FIC 
scheme with a time-step twice as long. 

In the general framework, when M. is not linear, these various schemes cannot be 
gathered in the unique formalism ((Zj) or ((HI). 

Addition of a second-order time-filter in the above definitions for 2-TL schemes is 
straightforward. Two main variants are usually considered, depending on whether the 
filtering is applied only to the time-averages of linear terms (as in SHB78), or also to 
the time-averages of non-linear terms, in (jH)). For instance, in a 2-TL SI scheme, the first 
variant of the filter consists in replacing £*(A'+(") + A'O) by /:*[(1 + k)A'+(") + (1-2k)A'° + 
kX~] in ((Hj), where k is a (small) positive parameter, and for the second variant, the same 
modification is also applied to the first RHS term of ((HJ. The scheme then becomes 
essentially a 3-TL scheme since information at level X~ is always used. However, the 
use of large values of k (e.g. k = 0.5, which eliminates the contribution) is known 
to deteriorate the solution through a spurious damping of transient perturbations (e.g. 
Hereil and Laprise, 1996). 

For 3-TL schemes, second-order time-filters are uneffective, and a first-order accurate 

time-decentering must be used. This consists in replacing C*(X'^^^^ + X^) by £*[(! + 

g-j^+(n) -|- (1 — 6)^:"^] in (|TTHl . where e is a (small) positive parameter. The scheme then 

ceases to be second-order accurate in time since the time-average is no longer centred in 

time. This results in a spurious damping of transient perturbations even for moderate 
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values of e (due to the weaker time-selectivity of the filter e compared to n). 



4 Conditions for space-continuous analyses 



4.1 Conditions on the upper and lower boundaries 



Space-continuous analyses are much easier to carry out when the equation system is 
defined in the whole unbounded atmosphere, because the expression of the normal modes 
of the system is more general. The following space-continuous analyses will restrict to this 
case (although this is not strictly required). For systems in which the vertical direction is 
represented (e.g. HPE and EE systems), this means that the upper and lower boundary 
conditions must not appear explicitly in the set of governing equations. However for 
systems cast in mass-based coordinates [such as HPE system in pressure-based coordinate 
(e.g. SHB78), and EE system in hydrostatic pressure-based coordinate (Laprise, 1992, 
L92 hereafter)], the upper and lower boundary conditions actually appear inside the set 
of equations through vertical integral operators with definite bounds at the boundaries of 
the vertical domain. When they are present, it is assumed that these integral operators 
can be eliminated, i.e. that C, C* can be transformed to "unbounded" operators by 
application of appropriate vertical linear differential operators to the prognostic equations 
which originally involve integral operators, in order that e.g. (jH)) rewrites as: 



d_ 

di 



\ 



IpXp 



( 



J 
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V 



PI 



IpCip 



IpjCpp 



\ 



J 



yXp J 



(11) 



where P is the number of prognostic variables of the unbounded system, {li, . . . ,lp) are 

linear vertical operators, and {liCu, . . . ,lpCpp) are linear spatial operators which no 

longer contain any reference to the upper and lower boundaries. The transformed system 
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obtained for L can then be written as: 

where / is the diagonal matrix (/i, . . . ,/p). A similar condition must be true for £* as 
well: it is assumed that applying the same operator / to £* leads to an operator IL* 
for which does not contain any reference to the upper and lower boundaries for 
(i, j) G (1, . . . , P). The first condition for the following analyses is: 

[CI]: There exists a linear operator / such as IL and /£* have no reference to the 
upper and lower boundaries. 

The system (|T2| is henceforth referred to as the "unbounded" system. 
4.2 Conditions on the stability of the X state 

The aim of these analyses is to determine in which conditions a stationary state X for 
IL will remain a stable equilibrium-state in the time-discretised context, provided it is 
a stable equilibrium-state in the time-continuous context. Hence the analyses will be 
restricted to stationary states X which are in stable equilibrium. Given the linear context 
used here, a physical transposition of this condition is: 

[C2]: For any perturbation X{t = 0) around X with a bounded energy-density, the 
time-evolution X{^) resulting from (fT2| must have a bounded energy-density. 

The condition is formulated with energy-density instead of total energy because the do- 
main is unbounded in space. The complex eigenmodes of the unbounded system fT^ are 
the complex functions of space X{v) which satisfy: 

iLXij) = XlX{r) (13) 
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where A G C and r denotes the spatial dependency. The time-evolution of the mode is 
bounded in time if A is a pure imaginary number, and Af(r) is then a "normal mode" in 
the usual terminology. Theoretical arguments out of the scope of this paper show that a 
mathematical transposition of [C2] writes: 

[C2']: For any complex eigenmode of /£, A G iH <^=^ '^(i") has a bounded 
energy-density. 

4.3 Conditions on the normal modes of the Hnear unbounded 
system 

The time-continuous normal modes of the original system ^ are the complex functions 
of space X{r) which have an oscillatory temporal evolution. Hence they satisfy CX{r) = 
iZJX{r). Similarly, the time-continuous normal modes of the linear unbounded system 
(fT2|) are the complex functions of space X{r) for which IX (r) has an oscillatory temporal 
evolution. Hence they satisfy l.CX{r) = iujl.X{r). 

Any normal mode of the original system is also necessarily a normal mode of the un- 
bounded system, with the same frequency U. For any normal mode of the unbounded 
system , we can choose the origin of space o in such a way that this mode writes: 

Xir) = X{o)f{r) = Xfir) (14) 

with X = {Xi, . . . , Xp) G C^, and / = (/i, . . . , fp) is a vector of space-dependent func- 
tions (the product Xf is understood "component by component"). For the indices i such 
as Xi uniformly vanishes, Xi = and /j(r) = 1 are assumed by convention. The function 
/ and the vector X are respectively termed the "structure" and the "polarisation vector" 
of the mode. The two following conditions are required for the proposed analyses: 

[C3]: For any normal mode X of the unbounded linear atmospheric system with a 

structure /(r), li.fiijc) must be proportional to /i(r): 
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Vz G (1, . . . , P) , hf,{r) = with G C*. (15) 

and: 

[C4]: For any normal mode X of the unbounded hnear atmospheric system with a 

structure /(r), Cij.fjijc) [resp. £*,-./j(r)] must be proportional to /j(r): 
V(z,j), kCij.fjir) =JIijfi{r) and liC*j.fj{r) = fi*jfi{r), with {flip G C. (16) 

As will be seen below, these latter two conditions have the important consequence that for 
any normal mode of the unbounded system, each individual time-discretised prognostic 
equation for A'j(r) becomes a scalar equation. This key ingredient makes the analysis 
straightforward for every member of the ICI class. 

4.4 Comments 

Since the set of normal modes for the unbounded system encompasses the set of normal 
mode of the bounded system, the transformation from the bounded system to the un- 
bounded system is not likely to "mask" some instabilities of the original system, unless the 
causes of the instability lie in the boundary conditions themselves. However, discretised- 
analyses allow to clarify this point by showing that the stability of the bounded and un- 
bounded systems are actually found to be similar in practice. Besides, discretised-analyses 
of the unbounded system are by nature impossible to perform, hence the continuous anal- 
ysis is the only way to estimate the intrinsic stability of the unbounded system and to 
demonstrate that instabilities, when they occur in a practical application, are not due to 
a weakness in the spatial discretisation or even to the boundary conditions, but actually 
to the time-discretised propagation of free modes inside the atmosphere. 

In spite of their apparently abstract and constraining form, conditions [C1]-[C4] are 

easy to verify with routine normal mode analysis techniques when examining a particular 
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concrete meteorological system. The condition [C2'] restricts the set of stationary states 
X around which the analysis is meaningful, and conditions [CI], [C3], [C4] restrict the 
spectrum of meteorological contexts accessible to the analysis, since they require a quali- 
tative similarity between /, C and C* operators. As stated in SHB78, analyses performed 
under this type of conditions "...grossly exaggerate the stability of the scheme..." since 
in more realistic meteorological contexts, the atmospheric and reference operators can be 
qualitatively much more different than imposed by this condition. 

5 Time-Discretised Space- Continuous Analysis 

The analysis examine the stability of the time-discretised system for perturbations which 
have a time-continuous normal mode structure. Hence we consider a given function / = 
(/i,...,/p) which is a normal mode structure for the time-continuous system, and we 
determine the normal modes of the time-discretised system which have the same structure 
/, by solving the equation: 

^{t=Ai)/(r) = Ai=(i=o)/(r) (17) 

where ^'(t^o) cind A are the unknowns, and Xi^t=At) is determined using the time-discretisation 
scheme ((Zj) or For schemes using three time levels (as Leap-Frog or extrapolating 2- 
TL) a similar relationship A:'(t=-At)/(r) = A~^A'(t=o)/(i") must be added. If for some 
solution, |A| > 1 (resp. < 1) the scheme is unstable (resp. damping) for this particular 
mode. The ratio Arg(A)/(zu7At) gives the relative phase-speed error of the scheme for 
this mode. 

The analysis is described here for a 2-TL discretisation ((Zj) , but the transformation to 

a 3-TL scheme as well as the addition of time-filters such as k or e are straightforward. 

In the following of this section, the notation X{r,t) is replaced by the usual superscript 

notation for time-discretised variables ^*(r) as in sectional 
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As a consequence of the discussion in section|Hland applying |T7l) , X'^^'^^ can be written 
as X~^^'^\r) = yu(A)A:'°(r), where /i(A) depends on the choice of the initial guess (e.g. 
/X = 1 for a 2-TL non-extrapolating scheme, /i = 2 — 1/Afora 2-TL extrapolating scheme, 
etc.). The original unbounded system (|T^ is thus time-discretised following 0: 



/A-+W(r) -/A'°(r) 
At 



+ ne 1, 



A^iter) 







For schemes with Nuer > 2, it is assumed that the intermediate states lX~^^'^\r) for 
i G (1, . . . , A'^iter — 1) have the same structure / as the one currently examined, which 



allows to define the polarisation vector X^^"'^ by: 



+ (n), 



(18) 



Hence, using (IT^ and (|TH|l . the space dependency /j(r) eliminates. For the generalized 
state-vector Z = (A^o, . . . , ^+(^''-)) of length (A^iter + 2) x P, the above 

system writes: 



v 



/i(A)/p -Ip Op 

Ml Ma M3 

: Op ••. 

Ml : 

-AJp Op ■ ■ ■ 



0, 



••. Op 

M2 M3 

Op Ip 



.Z = M.Z = 0. 



(19) 



where Ip (resp. Op) is the unit (resp. null) P-order matrix, and: 
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where 5ij is the j) Kronecker symbol. The possible values of A for the normal mode 
structure /(r) that we examine, are thus given by the roots of the following polynomial 
equation in A: 

Det(M) = (23) 

The dependencies to A are limited to the top- and bottom-left blocks. For a non- 
extrapolating 2-TL scheme the degree of the polynomial is P, and there are P phys- 
ical modes associated with this structure /(r). For time-schemes making use of three 
time-levels (i.e. 3-TL schemes, or extrapolating 2-TL schemes, or 2-TL schemes with a 
time-filter), the degree becomes 2P, and there are P additional computational modes. 
The growth-rate for any of these modes is given by the modulus of the corresponding 
complex root of (|^ . The growth rate of the time scheme for the considered structure / 
is then defined by the maximum value of the modulus of these P or 2P roots: 

r(/) = Max(|A,(/)|), zG(l,...,P) [or (1,...,2P)] (24) 

Analytical solution of (f23|) is not possible for large values of P, and a numerical solution 
is often needed. 

In this paper we will call "asymptotic growth-rate" the growth-rate for the matrix 
M in the limit of large time-steps At — » oo. The analysis of the asymptotic growth- 
rate is easier than for finite time-steps, since the matrix M of (f23|) degenerates to a 

matrix M' in which 5ij vanishes and At/2 eliminates. Another advantage of asymptotic 
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growth-rates is that they appear to be independent of the structure / in most of the cases 
examined below, thus simplifying considerably the interpretation of the results. When the 
asymptotic growth-rate is independent of the structure / the growth-rate of the scheme 
can be defined by the growth-rate obtained for this scheme with any structure /. 

When the growth rate for a given scheme is one (or less) for any mode of any normal 
mode structure /, the scheme is then said to be "unconditionally stable" (being under- 
stood "in At"). The criterion for unconditional stability obtained through Det(M') = 
is not only of academic interest since the considered time-schemes are actually used with 
large time-steps in NWP: the practice shows that a scheme which is not unconditionally 
stable in the simplified context of these analyses has few chance to be robust enough for 
use in real conditions. 



6 Simple examples: ID systems 
6.1 ID Shallow- water system 

The linearised ID shallow water system in an horizontal direction x can be classically 
written in terms of the wind u along x, and the geopotential (j): 

du d(j) 

m ^ ^^^^ 
Tt = -'^Tx 

This system is also valid for the external mode of an isothermal atmosphere in HPE and EE 

system, replacing by A{R^ /Cp)T (which is done in the following). The reference system 

is obtained by replacement of T by T*, and a "non-linearity" factor is defined through: 

a = (T - T*)/T*. Solution of ^ implies that if T < 0, X e iU =^ u = uexp{rx), 

which has not a bounded energy-density for rx — > +oo. Hence [C2'] requires T > (i.e. 

a > —1). The boundary conditions do not appear explicitly in the system, hence / can 
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be taken as the identity operator to satisfy [CI], and [C3]. In the notations of section 
m and El we have P = 2 , Xi = u and X2 = (p. The normal modes of the system write 
il){x) = i/jexp{ikx) with /c G R and i/j = {u,(f)). Conditions [CI] - [C4] are easily checked 
to be satisfied. 

For a 3-TL SI scheme, writes: 

/\2_1\2 1,2*2 

(^) =-^(A^ + l)(A^ + l + 2aA), (27) 

where c* = 2a/ {R/Cp)RT* . In the limit of long time-steps, the LHS term disappears, 
and the four roots of the RHS give the "asymptotic" numerical growth-rate for the two 
physical and two computational modes of the system. The two roots of the first factor 
have a neutral stability, while those of the second factor have a modulus equal to 1 if 
— 1 < a < 1. The criterion (on X, X*) for unconditional stability (in At) of the 3-TL 
SI scheme is thus : < T < 2T*. Some further algebraic manipulations from (f23|l with 
At = 00 show that this criterion remains unchanged when increasing A^'iter- 
For a 2-TL SI non-extrapolating (/i = 1) scheme, (f23|) becomes: 

[-^) =-^(A + l)(A + l + 2a), (28) 

and the criterion for unconditional stability becomes more constraining than for the 3- 
TL scheme: — 1 < a < (i.e. < T < T*). The asymptotic growth-rate of a 2-TL 
non-extrapolating ICI scheme with A'iter iterations is given by: 

r = Max (1, |2(-a)^"- - l|) (29) 

The domain for unconditional stability is thus — 1 < a < for odd values of Alitor, and 
— 1 < a < 1 for even values of A^'iter- 

6.2 ID vertical acoustic system in mass-based coordinate 

We consider a vertical ID compressible atmospheric column satisfying the conditions 

listed in sectional A regular mass-based coordinate a = (tt/tt^) is chosen following L92, 
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by making 1] = a, A{a) = and B{a) = cr in equation (31) of L92. The variable n denotes 
the hydrostatic pressure, and tt^, the surface hydrostatic pressure. The system is readily 
obtained from equations (36) - (45) of L92 by removing all horizontal dependencies. The 
surface hydrostatic-pressure does not evolve in time (see equation (45) in L92). 

The equations are linearised around a resting atmospheric-state X and a resting 
reference-state A"*, both satisfying the conditions of section [21 The temperatures T and 
T* are taken uniform, and we still define the "non-linearity" factor by: a = {T — T*)/T*. 
The pressure values p and p* are assumed to be equal to a common value ttq at the origin 
(cr = 1). Since X and X* are hydrostatically balanced, p = W = aiiQ and p* = n* = aiiQ at 
any level. The thermodynamics equation decouples and the linear system around X for 
the vertical velocity w , and the pressure deviation p' = p—p writes in standard notations: 

dw _ g_dj/_ 

dt ~ TToda 

dp' _ CpgTTo ^dw 
dt - C^Rt'' da 

The same derivation holds for C* and leads to an operator formally identical to the RHS 
of (f3fl|) - (f3T| . still acting on {w,p'), but with T replaced by T*. The solution of (fT3|l . 
implies that if T < 0, A G =^ w = w with r < — 1 or r > 0. For the mode with 
r < — 1, the energy-density is not bounded when a ^ 0. If T > 0, the structure of the 
normal modes of (f30|l - (f3T| is given by: 

w{a) = w a^'""-^'^^ = w h{a) (32) 
p\a) = p'a^^''+'/''> =p' Ma) (33) 

where is a real number, and they have a bounded energy-density. The condition [C2'] 
therefore requires T > (i.e. a > —1). Finally, [C4] is trivially checked to be satisfied. 



For a 3-TL SI scheme, (j23|) writes: 



(A^ + l) A^ + l-f^ , (34) 



2 At J 4i7*2 V J \ 1 + 
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where c* = ^J{Cp/C^)RT* and H* = RT*/g. Comparison of (HI and ^ shows that 

the stability of the ID vertical system for a is the same as the one of the previous 

shallow-water system for a' = — «/(l + a). Hence the criteria for unconditional stability 

directly follows from those of the previous case, by similarity arguments. The criterion 

for unconditional stability of the 3-TL SI scheme is a > (—1/2), i.e. T > (1/2)T*, and 

this criterion remains unchanged when increasing iVitcr- 

For a 2-TL SI non-extrapolating (/i = 1) scheme, becomes: 
'A-lV (i^2 + l/4)c*2 ,/ 2a 



and the criterion for unconditional stability (a > 0, i.e. T > T*), which is more con- 
straining than for the 3-TL SI scheme. For iterated 2-TL schemes, the criterion for 
unconditional stability is a > ( — 1/2) for even values of A'iter, and a > for odd values. 

6.3 ID vertical acoustic system in height-based coordinate 

In this example we show that the stability properties of the ID vertical system may depend 
on the coordinate. The framework is taken as in the previous example except that the 
vertical coordinate is the height z. The linearised system C (cf. e.g. Caya and laprise, 
1999) writes: 

T' 

+ 9^ (36) 
w (37) 

Is) »• P«) 

where q' = q — q, q = \n{p/po), q = —gz/RT, pq is a reference pressure, and p is the true 

pressure. The normal modes of C have the following form, for = [w, T', q'): 

1 



dw 
'dt 


dz 


dr 


RT d 


dt 


Cy dz 


dq' 


- 


'dt 


\RT 


, q = 


-gz/RT, po 



ip = ip exp 



lu H — 
■211 



(39) 
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where H = (RT/g). The reference system C* is defined in a similar way replacing T by 
T* and q hj q* = —gz/RT*. It should be noted that the structure of the normal modes 
of C* is not the same as for C since the characteristic height for C* is H* = {RT*/g) 
For a 2-TL SI non-extrapolating (/i = 1) scheme, (|23|) becomes: 

— — = — iz/ + ^ A + 1 + 2a iu(X + l)- — { A + 1 

\ At J 4\ 2hJ^ ^[^ ^ H*\l + aJ\ 1 + 2a J 

(40) 

where c* = ^JJC^^jQi^)RT* . In the height-coordinate framework, the asymptotic growth- 
rate depends on the structure since v appears in one of the factors which become dominant 
at large time-steps. The interpretation of the results is thus slightly complicated in com- 
parison to the case of a mass-based coordinate. For the most external structure (z/ = 0), 
the asymptotic growth-rate is given by the roots of: 

This polynomial consists in a combination of two factors similar to those obtained in 
the two previous examples (through a formal replacement of 2a by a for the second 
factor). As a consequence, the unconditional stability domains can be readily deduced 
from these previous examples: the external structure z/ = is unstable for any value 
a 7^ when At — > oo. The instability is thus much more severe than in the case of a 
mass-based vertical coordinate, for which a > was sufficient to ensure unconditional 
stability. Moreover, slightly shorter structures with vertical wavelengths of the order of 
(l/H) are found to be more unstable than the external one, and for these structures, the 
unconditional stability criterion (a = 0) remains unchanged when iViter is increased. Fig. 
[U shows the asymptotic growth-rates for the 2-TL SI (A'^itcr = 1) scheme and the the 2-TL 
ICI scheme with TViter = 2 for z/ = 0.0001 m~\ a structure for which the instability is 
close to its maximum. The severe instability of the 2-TL SI scheme is only alleviated but 
not suppressed by choosing iViter = 2. 
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For the 3-TL SI scheme, the external structure u = is unconditionally stable for 
—0.25 < a < 1, but slightly shorter structures as above are found unstable at large 
time-steps as soon as a 7^ (very short modes are stable however). Fig. [T] depicts 
the asymptotic growth-rates for two structures: the external structure u = 0, and a long 
structure u = 0.0001 m~^. The growth-rate of the long structure for a moderate time-step 
At = 30 s with a time-decentering e = 0.1 (as in Caya and Laprise, 1999) is also depicted: 
the practical instability becomes small in these conditions, and the 3-TL scheme cannot 
be positively rejected, especially considering the fact that dissipative processes could act 
in a way to stabilize the scheme. The practical impact of this predicted weak instability 
for NWP applications could be easily evaluated with a ^-coordinate model, using an 
experimental set-up similar to the one used here, and then progressively extending the 
set-up to approach real-case experimental conditions. 

6.4 Comments 

In the three simple examples examined above, the criterion for unconditional stability 
is seen to be more constraining in the 2-TL non-extrapolating SI scheme than in the 
3-TL SI scheme. The 2-TL extrapolating SI scheme is found to have similar domains 
of unconditional stability than its non-extrapolating counterpart (not shown). For mass- 
based coordinates, if both vertically propagating acoustic waves and external gravity waves 
are simultaneously allowed by a given equation system, the above analyses suggest that 
2-TL SI schemes are so constraining that there is no domain for unconditional stability. 
For height-based coordinates, the long vertically propagating acoustic waves are always 
unstable in the 2-TL SI scheme. 

This leads to suspect that, in opposition to 3-TL SI schemes, classical 2-TL SI schemes 

are not suitable for the EE system with any vertical coordinate. The HPE system with 
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2-TL SI scheme did not suffer from this problem since vertically propagating waves are not 
allowed in HPE (i.e. the ID column atmosphere is stationary). The intrinsic instability 
of the 2-TL SI scheme for EE system is confirmed in sectional for mass-based coordinates. 

When a second-order time-filter with parameter k is applied to the system examined 
in the first example for a 2-TL SI non-extrapolating scheme, (f28|) becomes: 



At 



{X + 1) + k[X-2 + j 



(A + 1 + 2a) + K ( A - 2 + ^ 



(42) 

and the criterion for unconditional stability becomes — 1 < a < 2k,. For the second exam- 
ple, the similarity argument shows that the criterion for unconditional stability becomes: 
a > —2k /{I + 2k). The domains of stability of 3-TL and 2-TL time-filtered ICI schemes 
for the two first examples are summarized in Table 1. 

The application of a time-filter thus allows to alleviate the stability constraints for 
2-TL SI schemes, and a non-vanishing domain for unconditional stability is recovered. 
The width of the unconditional stability domain increases with k. This is found to hold 
for the ID vertical system in height-based coordinates as well, which is consistent with 
the results of Semazzi et al. (1995) and Qian et al (1998): they succeeded to solve 
numerically the EE system at low-resolution with a 2-TL SI scheme, however, the use of 
a large value k = 0.5 was required to stabilize the model. As a consequence the forecasts 
suffered from a dramatic loss of energy with increasing forecast-range, and ceased to be 
of meteorological interest after 2-3 days. Moreover, at high resolutions (and consequently 
steep orography) the use of a time-filter k is found experimentally to be an insufficient 
solution for eliminating the intrinsic instability of the scheme (not shown). 

If a high level of accuracy is desired for the EE system with a 2-TL classical SI time- 
discretization and high resolution, a more robust scheme (e.g. with a larger value of 

A^iter) must be used. The above analyses show that the unconditional stability domain is 
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dramatically reduced for odd values of iVjter, hence ICI schemes with even values of iViter 
are preferable for solving the EE system with a 2-TL scheme. 

The ID vertical system in mass-based coordinates has been found to be more stable 
than its counterpart in height-based coordinates in a general way. For mass-coordinates, 
the 3-TL SI and the 2-TL ICI schemes with even values of A'^iter have an extended domain 
of unconditional stability, whilst for height-based coordinates, they are unstable as soon 
as a 7^ 0. For 3-TL SI schemes, the necessity to have recourse to a first-order time- 
decentering e > to overcome this instability is a significant drawback since it results in a 
spurious damping of transient phenomena, similarly to k but in an even less selective way, 
as mentioned above. We think these differences give a substantial theoretical advantage 
to mass-based coordinates for solving the EE system with classical SI or ICI schemes. 

7 Analysis of the EE system for isothermal atmospheres 
in mass-coordinate 

The analysis of the isothermal HPE system for ICI schemes does not substantially modify 
the general conclusions drawn for the shallow- water case (not shown), hence the case of 
the EE system is directly examined. 

In this section, the EE system is cast in the pure unstretched terrain-following co- 
ordinate a which can be classically derived from the hydrostatic-pressure coordinate tt 
of L92, through a = (vt/tts) G [0,1], where vr^ is the hydrostatic surface-pressure. The 
nonhydrostatic prognostic variables are the non-dimensionalised nonhydrostatic pressure 
departure V = {p — n)/n (where p is the true pressure), and the vertical divergence d 
which writes in cr coordinate: 



d 



RT 



9 



dw 
drj 



(43) 
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The adiabatic system writes: 



r/V ET / cTP\ 

^ = -RTVq--^^VV-[l + V + a^]V(t) (44) 



dt ^ {l + V) V da J 



dt RT \ da J \ da 

9{l + V) 



d(VV - D, 



3) 



RT 

dT RT 



Vwia-^ 



(45) 



dt a"^^ ^''^ 



9g 

at 



- / (yV + VVq)da' 
Jo 



(48) 



where: 



Ds = VV + d + ^^^-^ V0 Ya^^ (49) 



i2T da J 

.-.f(T^)^ 

- = VVg - - / (VV + VVq)da', (51) 
71" cr Jo 



V is the horizontal wind, and V is the horizontal derivative operator. The domain is 
restricted to a vertical plane along {x, a) directions for clarity. The system is linearized 
around a resting isothermal and hydrostatically-balanced state X: 



-RQV^T + RT{g - J) V^P - RTV^q (52) 



dD 
'dt 



dt RT \ da J \ da 

dT RT 



dt a 



(D + d) (54) 



V 



dq 
di 



-AfD, (56) 
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where the vertical integral operators Q, S and ^f are defined by: 

gX = [ {X/a')da' (57) 



SX = (1/ct) / Xda' (58) 
Jo 

MX = [ Xda' (59) 
Jo 

The C* operator is similar to the RHS of this system, simply replacing T by T*. 
7.1 Verification of conditions [CI] - [C4] 

The linear operator li = a{d/da) is applied to (|^ . and U = [I + a{d/da)] to (fHK|l . The 
q equation (f56|l decouples and we obtain a linear unbounded system, in which (f52|) and 
(|55|) are replaced by: 

- RV'T-Ht[.^^t)v^V (60, 

Hence we have P = 4, A" = {D, d, T, P). Using of the same operators (/i, 1^), the reference 
operator is also made free of any reference to the upper and lower boundary conditions, 
which shows that the condition [CI] is satisfied. Solution of (fT3|) shows that [C2] requires 
T > (i.e. a > —1). The normal modes of the system are then: 

^Ij{x, a)=ip exp{ikx) a^""'^/^^ (62) 

where {k, z/) G R and ip represents D, d, T or V. In this particular case, the /i, • • • , /4 
functions are all identical. The verification of [C3], [C4] proceeds easily, as in previous 
sections. 



7.2 Results 

As a first illustration of the results, the 
schemes are shown in Fig. El as function 



growth-rates of 2-TL non-extrapolating ICI 

of a = (T — T*)/T* with a moderate time- 
26 



step At = 20 s, for three particular mode structures: 

(i) an external mode {k = 0.0005 m~^, z/ = 0) 

(ii) a vertical very internal mode (A; = 0, z/ = 100) 

(ii) an intermediate slantwise mode {k = 0.0005 m^^, u = 3) 

The suspicions raised in the simple ID examples are confirmed: the internal vertically- 
propagating mode is unstable for a < while the external gravity mode is unstable for 
a > 0. Moreover, intermediate, slantwise-propagating modes are unstable in the whole 
domain, and the acoustic external mode (Lamb-wave) appears to be unstable for a < 
as well. The domain of stability vanishes, which confirms that 2-TL SI scheme is not 
relevant for solving the EE system. The effect of introducing a time-filter to save the 
situation is discussed below. 

The asymptotic growth-rates resulting from the EE system for At = oo are now 
examined. Similarly to most previous cases, they are independent of the geometry {k, 
u) of the mode. Fig. [HI shows the asymptotic growth-rates as a function of a for 2-TL 
non-extrapolating ICI schemes with Alitor = (1,2,3,4). As stated above, the SI scheme 
(-^itcr = 1) is unstable for any value of a. For even values of A'itci., the scheme has an 
"optimal" domain of unconditional stability —1/2 < a < 1, while for odd values, the 
scheme is unstable for all values of a. 

For 3-TL ICI schemes the domain of unconditional stability is— l/2<a<l indepen- 
dently of the values of A'^itor and k. The curves (not shown) are similar to those obtained 
for even values of A'^iter for 2-TL ICI schemes. 

The impact of applying a time-filter k = 0.1 to 2-TL ICI schemes is depicted on Fig 

m for the second variant (the first variant behave qualitatively in the same way) . The 

global impact is to "lower" the curves of the asymptotic growth-rates, and consequently, 

to increase the width of the unconditional-stability domain. However, large values of k 
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(e.g. K ^ 0.5) are required in order to obtain a wide stability domain, especially for the 
2-TL SI scheme, and this strategy is known to be irrelevant for NWP purpose. 

Finally, it is worth noting that the results obtained for the EE system are fully com- 
patible with the conclusions that can be drawn from the intersection of the domains of 
unconditional stability in Table 1 for the two simple frameworks examined above in mass- 
coordinate. The ability of these very refined frameworks to capture the essence of the 
behaviour of the time-discretised EE system in the limit of long time-steps makes them 
very useful tools to fully understand the underlying causes of its stability or instability. 

8 Conclusion 

A general method for investigating the stability of the ICI class of time-discretisations on 
canonical problems with various space-continuous equation systems has been presented. 
These ICI schemes are based on a separation of evolution terms between a simple linear 
operator and "non-linear" residuals. The method has been validated by confirming earlier 
results, then the application to new frameworks (equation systems or time-discretisation 
schemes) allowed to extend these results. The main conclusions drawn from this study 
are: 

(i) Even on very simple (ID) examples, the stability properties of time-discretisations 
for a given equation system are very dependent on fundamental choices (as e.g. the 
choice of the vertical coordinate). Hence, the import of conclusions drawn from a 
given analysis must be carefully limited to the examined framework. 

(ii) For the EE system, height-based coordinates have a theoretical disavantage com- 
pared to mass-based coordinates since they exhibit an intrinsic instability for (long) 
vertically propagating waves. 
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(iii) The 2-TL SI scheme is found not to be appropriate for the EE set of equations, what- 
ever coordinate is employed (using a time-filter results in an unacceptable degrada- 
tion of the solution). 

(iv) For the EE system, the 2-TL scheme with Alitor = 2 brings a dramatic increase of 
the stability compared to the 2-TL SI scheme (A'itcr = !)• This statement holds for 

even values of A'ltcr, while odd values leads to a significantly weaker stability. 

(v) As a consequence of the latter point, the 2-TL ICI scheme with A'^itcr = 2 seems 

worth to be considered for the EE system. 

However, as mentioned in SHB78, the stability inferred from this type of analyses is 
overestimated, and fiows in which the non-linearity comes from other sources than the 
discrepancy between the atmospheric and reference temperature profiles could reveal new 
instabilities in practice. For instance, in spite of its apparent "optimal" stability in the 
simplified context of this paper, the 3-TL SI scheme has proved to be not stable enough 
for solving numerically the EE system in realistic highly non-linear conditions at high 
resolutions, due to other terms treated explicitly. This point clearly demonstrates the 
limitations of this type of academic exercise. 

Nevertheless, in spite of its necessary limitations, this study can serve to distinguish 
schemes which are definitely not relevant for practical use from the others, and give a first 
theoretical justification for those which are worth considering. In agreement with Cote et 
al (1998) and Cullen (2000) we think that ICI schemes with Nit^j. > 2 are among the most 
approriate alternatives for integrating the EE system in highly non-linear conditions at 
fine-scale, including from the point of view of efficiency. 
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List of Tables 

Table 1: Domains of unconditional stability for time-filtered schemes for the two first ID 
examples. 
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Shallow-water 


ID vertical (mass) 


3-TL ICI 


-1 < a < 1 


-1/2 < a 


2-TL ICI (A^iter even) 


-!<«<! 


-1/2 < a 


2-TL ICI (A^iter odd) 


, / AT 

-1 <a < (2k) ^/^'t- 


_2ft;(l/Mtor) ^ 
1 + 2ft:(l/Mte.) ^ " 



Table 1: Domains of unconditional stability for time- filtered schemes for the two first 
ID examples. 
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List of Figures 

Fig. 1: Asymptotic growth-rates T for ID vertical system in z coordinate as a function 
of the nonlinearity parameter a. thin line: long mode {u = 0.0001 m~^) with 2-TL SI 
scheme; thick line: long mode with 2-TL ICI scheme A'ltcr = 2; dotted line: external 
mode {u = 0)with 3-TL SI scheme; dashed line: long mode with 3-TL SI scheme. Circles: 
practical growth-rate of 3-TL SI scheme for the long mode with At = 30 s and e = 0.1. 

Fig. 2: Growth-rate F with At = 20 s for the EE system with a 2-TL SI scheme as 
a function of the nonlinearity parameter a. solid line: external mode (i); dashed line: 
slantwise mode (ii); dot-dashed line: internal mode (iii). The left-part of the solid line 
represents an acoustic external mode. 

Fig. 3: Asymptotic growth-rates F for EE system with 2-TL ICI scheme as a function 
of the nonlinearity parameter a. solid line: TViter = 1; dashed line: TVjter = 2; dot-dashed 
line: Alitor = 3; dotted line: Alitor = 4. 

Fig. 4: Same as Fig 3, but with a time-filter k = 0.1. solid line: Ajtcr = 1; dashed 
line: A^iter = 2; dot-dashed line: Ajter = 3; dotted line: Ajter = 4. 
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Same as Fm 3. but with a time-filter k. = 0.1. solid line: M.™ = 1; dashed line: A^^.„, = 2; dot-da,f 
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Figure 1: Asymptotic growth-rates F for ID vertical system in z coordinate as a function 
of the nonlinearity parameter a. thin line: long mode {y = 0.0001 m~^) with 2-TL SI 
scheme; thick line: long mode with 2-TL ICI scheme Alitor = 2; dotted line: external 
mode {u = 0)with 3-TL SI scheme; dashed line: long mode with 3-TL SI scheme. Circles: 
practical growth-rate of the 3-TL SI scheme for the long mode with At = 30 s and e = 0.1. 
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Figure 2: Growth-rate T with At = 20 s for the EE system with a 2-TL SI scheme as 
a function of the nonlinearity parameter a. solid line: external mode (i); dashed line: 
slantwise mode (ii); dot-dashed line: internal mode (iii). The left-part of the solid line 
represents an acoustic external mode. 
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Figure 3: Asymptotic growth-rates T for EE system with 2-TL ICI scheme as a function 
of the nonlinearity parameter a. solid line: A'lter = 1; dashed line: Alitor = 2; dot-dashed 
line: iVjter = 3; dotted line: Alitor = 4. 
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Figure 4: Same as Fig 3, but with a time-filter k = 0.1. solid line: iViter = 1; dashed line 
-^iter = 2; dot-dashed line: iViter = 3; dotted line: iViter = 4. 
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